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ABSTRACT 



Magnetised turbulence is ubiquitous in astrophysical systems, where it no- 
toriously spans a broad range of spatial scales. Phenomenological theories of 
MHD turbulence describe the self-similar dynamics of turbulent fluctuations in 
the inertial range of scales. Numerical simulations serve to guide and test these 
theories. However, the computational power that is currently available restricts 
the simulations to Reynolds numbers that are significantly smaller than those in 
astrophysical settings. In order to increase computational efficiency and, there- 
fore, probe a larger range of scales, one often takes into account the fundamental 
anisotropy of field-guided MHD turbulence, with gradients being much slower 
in the field-parallel direction. The simulations are then optimised by employing 
the reduced MHD equations and relaxing the field-parallel numerical resolution. 
In this work we explore a different possibility. We propose that there exist cer- 
tain quantities that are remarkably stable with respect to the Reynolds number. 
As an illustration, we study the alignment angle between the magnetic and ve- 
locity fiuctuations in MHD turbulence, measured as the ratio of two specially 
constructed structure functions. We find that the scaling of this ratio can be 
extended surprisingly well into the regime of relatively low Reynolds number. 
However, the extended scaling becomes easily spoiled when the dissipation range 
in the simulations is under-resolved. Thus, taking the numerical optimisation 
methods too far can lead to spurious numerical effects and erroneous represen- 
tation of the physics of MHD turbulence, which in turn can affect our ability 
to correctly identify the physical mechanisms that are operating astrophysical 
systems. 
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Subject headings: magnetic fields — magnetohydrodynamics — turbulence 



Introduction 



Magnetised turbulence pervades the universe. It is likely to play an important role 
in the transport of energy, momentum and charged particles in a diverse range of astro- 
physical plasmas. It is studied with regards to its influence on the generation of mag- 
netic fields in stellar and planetary interiors, small-scale structure and heating of stel- 
lar winds, the transport of angular momentum in accretion discs, gravitational collapse 



and interstellar scintillation (e.g.. 


Bis 


camp 


2OO3I: 


Kulsrud 


2004 




McKee & Ostriker 


2007 




Gold 


stein, Roberts & Matthaeus 


1995; 


Brandenburg & Nordlund 


2011; 


Schekochihin & Cowlev 


2007 


). The effects of magnetised turbulence need to be taken into account when analysing 



astrophysical observations and also when modelling astrophysical processes. 

The simplest theoretical framework that describes magnetised plasma turbulence is that 
of incompressible magnetohydrodynamics (MHD), 

d 



dt 



T ■ V z± + (z^ ■ V) z± = -VP + z/V^z^ + f ^ 



V ■ z± = 0, (2) 

where the Elsasser variables are defined as z^ = v ± b, v is the fluctuating plasma velocity, 
b is the fluctuating magnetic field normalized by y^AjipQ, = Bo/v^vrp^ is the Alfven 
velocity based upon the uniform background mag netic field Bq, P = {p/po + p is the 

plasma pressure, po is the background plasma density, represents forces that drive the 
turbulence at large scales and for simplicity we have taken the case in which the fluid viscosity 
u is equal to the magnetic resistivity. Energy is transferred t o smaller scales by the nonlinear 
interactions of oppositely propagating Alfven wavepackets (IKraichnanl Il965[ ) . This can be 
inferred directly from equation ([T]) by noting that in the absence of forcing and dissipation, 
if z'^{'x,t) = then any function z''=(x, t) = F''=(x ± V^t) is an exact nonlinear solution 
that propagates parallel and anti-parallel to Bq with the Alfven speed. The efficiency of 
the nonlinear interactions splits MHD turbulence into two regimes. The regime in which 
the linear terms dominate over the nonlinear terms is known as 'weak' MHD turbulence, 
otherwise the turbulence is 'strong'. In fact, it has been demonstrated both analytically and 
numerically that the MHD energy cascade occurs predominantly in the plane perpendicular 
to the guiding magnetic field. This ensures that even if the turbulence is weak at large scales 
it encounters the strong regime as the cascade proceeds to smaller scales. MHD turbulence 
in astrophysical systems is therefore typically strong. 
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For strong MHD turbulence, iGoldreich fc Sridharl (jl995[ ) argued that the hnear and 
nonhnear terms in equations ([1]) should be ap proximately balanced at all scales, known as 
the critical balance condition. Consequently, IGoldreich fc Sridharl ( 19951 ) postulated that 
the wave packets get progressively elongated in the direction of the guide field as their 
scale decreases (with the field-parallel lengthscale / and field-perpendicular scale A related 
by / oc A^/^) and that the field-perpendicular energy spectrum takes the Kolmogorov form 
EGsik±) oc k] 



-5/3 



Recent high resolution direct numerical simulations with a strong guide field {va > 



5Vr 



perpendicular energy spectrum appears to be closer to E(k \ ) oc k i^^"^ (e.g., 


Maron & Goldreich 


2001: Miiller & GraDDin 2005: Mason et al. 


2008: 


Perez & Boldvrev 2009 


2010 


). A resolu- 


tion to this contradiction was proposed in 


Boldvre\ 


( 


2006f). Therein it was sueeested that in 



addition to the elongation of the eddies in the direction of the guiding field, the fiuctuating 
velocity and magnetic fields at a scale A ~ l/k± are aligned within a small scale-dependent 
angle in the field perpendicular plane, 9x oc A^/^. In this model the wavepackets are three- 
dimensionally anisotropic. Scale-dependent dynamic alignment reduces the strength of the 
nonlinear interactions and leads to the field-perpendicular energy spectrum E{k±) oc k^^^'^. 

Although the two spectral exponents —5/3 and —3/2 are close together in numerical 
value, the physics of the energy cascade in each model is different. The difference between 
the two exponents is especially important for inferring the behaviour of processes in astro- 
physical systems with extended inertial intervals. For example, the two exponents can lead 
to noticeably di fferent predictions for the rate of turbul e nt heating in coronal holes and the 



solar wind (e.g., IChandran et al.ll2010l : IChandranI l2010t iPodesta fc Borovskyll2010l ). Thus 



there is much interest in accurately determining the spectral slope from numerical simula- 
tions. Unfortunately, the Reynolds numbers that are currently accessible by most direct 
numerical simulations do not exceed a few thousand, which complicates the precise identifi- 
cation of scaling exponents. Techniques for careful optimisation of the numerical setup and 
alternative ways of differentiating between the competing theories are therefore much sought 
after. 

Maximising the extent of the inertial range is often achieved by implementing phys- 
ically motivated simplifying assumptions. For example, since the turbulent cascade pro- 
ceeds predominantly in the field-perpendicular plane it is thought that the shear-Alfven 
waves control the dynami cs while the pseudo-Alfven waves play a passive role (see, e.g., 
Maron fc GoldreichI ( 1200 ll )). If one neglects the pseudo-Alfven waves (i.e. removes the 
fiuctuations parallel to the strong guide field) one obtains a system that is equivalent 
to the reduced MHD system (RMHD) that was originally derived in the context of fu- 
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sion devices by iKadomtsev fc Pogutsd (119741 ) and [Strauss Jl976h ( see also iBiskampI ( 120031 ) ; 
Oughton. Dmitruk fc MatthaeusI (120041 ) ). Incompressibility then enables the system to be 
further reduced to a set of two scalar equations for the Elsasser potentials, resulting in a 
saving of approximately a factor of two in computational costs. Further computational sav- 
ings can be made by making use of the fact that the wavepackets are elongated. Hence 
variations in the field-parallel direction are slower than in the field-perpendicular plane 
and a reduction in the field-parallel resolution would seem possible. Indeed, this is widely 
used as an optimisati on tool in numerical simulations of the iner t ial range of field-guided 



MHD turbulence (e.g..lMaron fc Goldreichll200ll : iMason et al.ll2008l : iGrappin fc Mullerll2010 
Perez fc BoldvrevI l201of )~ The accumulated computational savings can then be re-invested 



in reaching larger Reynolds numbers for the field-perpendicular dynamics. 

Additionally, it is advantageous to seek other ways of probing the universal scaling of 
MHD turbulence. In this work we point out a rather powerful method, which is based on the 
fact that there may exist certain quantities in MHD turbulence that exhibit very good scaling 
laws even for turbulence with relatively low Reynolds numbers. The situation here is reminis- 
cent of the well kn own phenomenon of extended self-similarity in hydrodynamic turbulence 
( iBenzi et al.lll993l ). We propose that one such "stable" object is the alignment angle between 
the velocity and magnetic fluctuations, which we measure as the ratio of two specially con- 
structed structure functions. This ratio has been recently measured i n numerical simulations 
i n an a ttempt to differenti ate among various theoretical predictio ns (IBeresnvak fc Lazarian 
( 120061 ): iMason et al.l ( 120061 )). Also, it has recently been shown by iPodesta et al.l ( 120091 ) that 
the same measurement is accessible through direct observations of solar wind turbulence. 
Scale- dependent alignment therefore has practical value: its measurement may provide an 
additional way of extracting information about the physics of the turbulent cascade from 
astrophysical observations. In the present work we conduct a series of numerical simulations 
with varying resolutions and Reynolds numbers. We find that as long as the simulations 
are well resolved, the alignment angle exhibits a universal scaling behavior that is virtually 
independent of the Reynolds number of the turbulence. Moreover, we find that the length of 
scaling range for this quantity extends to the smallest resolved scale, independently of the 
Reynolds number. This means that although the dissipation spoils the power-law scaling 
behaviour of each of the structure functions, the dissipation effects cancel when the ratio of 
the two functions is computed and the universal inertial-range scaling extends deep in the 
dissipation region. 

The described method allows the inference of valuable scaling laws from numerical 
simulations, experiments, or observations of MHD turbulence with limited Reynolds number. 
However, one can ask how well the extended-scaling method can be combined with the 
previously mentioned optimisation methods relying on the reduced MHD equations and 



- 5 - 



a decreased parallel resolution. We check that reduced MHD does not alter the result. 
However, when the dissipation region becomes under-resolved (as can happen, for example, 
when the field-parallel resolution is decreased) , the extended scaling of the alignment angle 
deteriorates significantly. Thus the optimisation technique that works well for viewing the 
inertial range of the energy spectra should not be used in conjunction with the extended- 
scaling measurements that probe deep into the dissipation region. The remainder of this 
paper will report the findings of a series of numerical measurements of the alignment angle 
in simulations with different Reynolds numbers and different field-parallel resolutions in both 
the MHD and RMHD regimes. The aim is to address the need to find an optimal numerical 
setting for studying strong MHD turbulence and to raise caution with regards to the effects 
that implementing simplifying assumptions in the numerics can have on the solution and its 
physical interpretation. 



2. Numerical results 

We simulate driven incompressible magnetohydrodynamic turbulence in the presence of 
a strong uniform background magnetic field, Bq = BqBz. The MHD code solves equations 
f lT]l2|) on a periodic, rectangular domain with aspect ratio x Ly, where the subscripts 
denote the directions perpendicular and parallel to Bq and we take L±_ = 27t. A fully 
dealiased 3D pseudospectral algorithm is used to perform the spatial discretisation on a 
grid with a resolution of Nj_ x N\\ mesh points. The RMHD code solves the reduced MHD 



counterpart to equations ([HE]) in which — (z^,z^, 0) fsee iPerez fc BoldvrevI (120081 )1 



The domain is elongated in the direction of the guide field in order to accommodate the 
elongated wavepackets and to enable us to drive the turbulence in the strong regime wh ile 



maintaining an inertial range that is as extended as possible (see lPerez fc BoldyrevI fl2010l )). 
The random forces are applied in Fourier space at wavenumbers 2tt/L± < fcj, < 2(27r/Lj_), 
(27r/L||) < k\\ < (27r/L||)r2||, where we shall take n\\ = 1 or n\\ = 2. The forces have no 
component along z and are solenoidal in the xy-plane. All of the Fourier coefficients outside 
the above range of wavenumbers are zero and inside that range are Gaussian random numbers 
with amplitudes chosen so that Vrms ~ 1- The individual random values are refreshed 
independently on average every r = nrLj_/27iVrms, i-e. the force is updated approximately 
l/ur times per turnover of the large-scale eddies. The variances cx^ = (If"*^!^) control the 
average rates of energy injection into the z~^ and z~ fields. The results reported in this paper 
are for the balanced case (T+ ~ cr_. In all of the simulations performed in this work we will 
set the background magnetic field i?o = 5 in velocity rms units. Time is normalised to the 
large scale eddy turnover time tq = L±/27iVrms- The field-perpendicular Reynolds number 



- 6 - 



is defined as Re±_ = Vrms{L±/27r)/u. 

The system is evolved until a stationary state is reached, which is confirmed by observing 
the time evolution of the total energy of the fluctuations, and the data are then sampled in 
intervals of the order of the eddy turnover time. All results presented correspond to averages 
over approximately 30 samples. We conduct a number of MHD and RMHD simulations with 
different resolutions, Reynolds numbers and field-parallel box sizes. The parameters for each 
of the simulations are shown in Table [TJ 

For each simulation we calculate the scale-dependent alignment angle between the shear- 
Alfven velocity and magnetic field fiuctuations. We therefore define velocity and magnetic 
differences as = v(x + r) — v(x) and 5b,. = b(x + r) — b(x), where r is a point-separation 
vector in the plane perpendicular to Bq. In the MHD case the pseudo-Alfven fluctuations 
are removed by subtracting the component that is parallel to the local guide fleld, i.e. we 
construct ^v^ = ^v^ — {Svr ■ n)n (and similarly for 6hr) where n = B(x)/|B(x)|. In the 
RMHD case fluctuations parallel to Bq are not permitted and hence the projection is not 
necessary. We then measure the ratio of the second order structure functions 

(|5Vr||5br|) 

where the average is taken over different positions of the point x in a given fleld-perpendicular 
plane, over all such planes in the data cube, and then over all data cubes. By deflnition of 
the cross product 

^ ^ sin(^,) ^ ^„ (4) 
where Or is the angle between and 6hr and the last ap proximation is valid for small 



angles. We recall that the theoretical prediction is 9r oc r^/^ f BoldyrevI I2OO6I ) . 



Figure [T] illustrates the ratio (jl]) as a function of the separation r = |r| for two MHD 
simulations (Ml and M2) corresponding to a doubling of the resolution from 256^ to 512'^ 
mesh points with the Reynolds number increased from Re± = 800 to 2200. Excellent agree- 
ment with the theoretical prediction 9 oc r^/^ is seen in both cases. As the resolution and 
Reynolds number increase, the scale-dependence of the alignment angle persists to smaller 
scales. Indeed, we believe that the point at which the alignment saturates can be identi- 
fled as the dealiasing scale, kd = N/3 = 85, 170 corresponding in conflguration space to 
~ l/2kii ^ 0.006,0.003 for the 256^,512^ simulations, respectively. This is verifled in 
Figure [2] that shows that alignment is largely insensitive to the Reynolds number (provided 
that the system is turbulent) and Figure [3] that shows that the saturation point decreases 
by a factor of approximately 2 as the resolution doubles at flxed Reynolds number. Thus as 
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computational power increases, allowing higher resolution simulations to be conducted, we 
expect to find that scale-dependent alignment persists to smaller and smaller scales. 

The fact that even in the lower Reynolds number cases scale-dependent alignment is 
clearly seen over quite a wide range of scales is particularly interesting, as in those cases 
only a very short inertial range can be identified in the field-perpendic ular energy spectrum , 



sp< 



making the identification of spectral exponents difficult (see Figure 1 in lMason et al.l ( 120061 )). 
In the larger Re cases, we can estimate the inertial range of scales in configuration space to be 
the range of r ~ l/2k over which the energy spectrum displays a power law dep endence. The 



f ield-p erpendicular energy spectrum for the Case M2 is shown in Figure 1 in iMason et al. 



( 120081 ). with the inertial range corresponding to approximately 4 < A; < 20, i.e. 0.025 ^ r < 
0.125. Comparison with Figure [T] shows that a significant fraction of the region over which 
the scaling 9^ oc r^/^ is observed corresponds to the dissipative region, i.e. that ratios of 
structure functions appear to probe deeper than the inertial range that is suggested by the 
energy spectra. 

We now consider the effect on the alignment ratio of decreasing the field-parallel resolu- 
tion. Figure m shows the results from three MHD simulations (M3, M4 & M5) for which the 
field-parallel resolution decreases by a factor of two, twice. As the resolution decreases the 
extent of the self-similar region diminishes and the scale-dependence of the alignment angle 
becomes shallower. If one were to calculate the slope for the lowest field-parallel resolution 
case (M5) one would find a scale-dependence that is shallower than the predicted power law 
exponent of 1/4. This may lead one to conclude (incorrectly) that scale-dependent alignment 
is not a universal phenomenon in MHD turbulence. However, the effect is obviously a result 
of the poor resolution rather than being an attribute of the alignment mechanism itself. 

Finally, we mention that for the three cases illustrated in Figure HJ the field-perpendicular 
energy spectra (not shown) display no appreciable difference. Since the Reynolds number is 
moderate the inertial range in /c-space is quite short. However, when the spectra are compen- 
sated with k^^'^ and k^^^ the former results in a better fit in all cases. This happens for two 
reasons. First, the stronger deviation from the alignment scaling 9 oc r^/^ occurs deeper in 
the dissipation region, that is, further from the inertial interval where the energy spectrum is 
measured. Second, according to the relationship between the scaling of the alignment angle 
and the energy spectrum, a noticeable change in the scaling of the alignment angle leads to 
a relatively small change in the scaling of the field-perpendicular energy spectrum. 
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Fig. 1. — The alignment angle for two MHD simulations (Ml and M2) corresponding to 
doubling the resolution and increasing the Reynolds number. In all the figures, the straight 
line has a slope of 1/4. 

3. Conclusion 

There are two main conclusions that can be drawn from our results. The first is that 
the measurement of the alignment angle, which is composed of the ratio of two structure 
functions, appears to display a self-similar region of significant extent, even in the moderate 
Reynolds number case which requires only a moderate resolution. We have checked that 
plotting the numerator and the denominator of the alignment ratio separately as functions 
of the increment r displays only a very limited self-similar region, from which scaling laws 
cannot be determined. A clear scaling behaviour is also not found whe n one plots the nu - 



merator versus the denominator as is the case in extended self-similarity (IBenzi et al.lll993l ). 
The result is interesting in its own right. It also has important practical value as it allows us 
to differentiate effectively between competing phenomenological theories through numerical 
simulations conducted in much less extreme parameter regimes than would otherwise be 
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Table 1. Simulation Parameters 



Case 


Regime 




iV|| 






rir 


Ml 


MHD 


256 


256 


5 


800 


1 


M2 


MHD 


512 


512 


5 


2200 


1 


M3 


MHD 


512 


512 


5 


2200 


0.1 


M4 


MHD 


512 


512 


10 


2200 


0.1 


M5 


MHD 


512 


256 


10 


2200 


0.1 


Rl 


RMHD 


512 


512 


6 


960 


0.1 


R2 


RMHD 


512 


512 


6 


1800 


0.1 


R3 


RMHD 


256 


256 


6 


960 


0.1 




Fig. 2. — The alignment angle for two RMHD simulations (Rl & R2) with the same resolution 
but different Reynolds number. 
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Fig. 3. — The alignment angle for two RMHD simulations (R2 & R3) with the same Reynolds 
number but different resolution. 
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Fig. 4. — The alignment angle for three MHD simulations (M3, M4 & M5) with different 
field-parallel resolutions 6\\ = {L\\/2'k) /N\\ and Re = 2200. 

necessary. The result could be especially useful if it extends to ratio s of structure functions 



for which an exact relation, such as the iPolitano fc Pouquetl (119981 ) relations, is known for 



one part, as it would then allow the inference of the scaling of the other structure func- 
tion. Reaching a consensus on the theoretical description of magnetised fluctuations in the 
idealised incompressible MHD system represents the first step towards the ultimate goal of 
building a theoretical foundation for astrophysical turbulence. 

The second main result that can be drawn from our work is that the measurement of the 
alignment angle appears to probe deep into the dissipation region and hence it is necessary 
to adequately resolve the small scale physics. As the field-parallel resolution is decreased, 
numerical errors contaminate the physics of the dissipative range and affect measurement 
of the alignment angle. As the decrease in resolution is taken to the extreme, the errors 
propagate to larger scales and may ultimately spoil an inertial range of limited extent. We 
propose that similar contamination effects should also arise through any mechanism that 
has detrimental effects on the dissipative physics. Mechanisms could include pushing the 
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Reynolds number to the extreme or using hyper diffusive effects. For exam ple, our results 



may provide an explanation for the numerical findings by iBeresnyakI (120111 ) who noticed a 
flattening of the alignment angle in simulations of MHD turbulence with a reduced parallel 
resolution and strong hyperdiffusivity. 

We also point out that the resul t recalls the phen omenon of extended self-similarity 



in isotropic hydrodynamic turbulence (IBenzi et al.lll993l ). which refers to the extended self- 



similar region that is found when one plots one structure function versus another, rather than 
as a function of the increment. Our finding is fundamentally different however, in the sense 
that the self-similar region only becomes apparent when one plots ratios of structure functions 
versus the increment, rather than structure functions versus other structure functions. Our 
result appears to be due to a non-universal features in the amplitudes of the functions, rather 
than their arguments, cancelling when the ratios are plotted. Whether such a property holds 
for other structure functions in MHD turbulence is an open and intriguing question. This is 
a subject for our future work. 
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